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ABSTRACT 

A computer technique for obtaining near-field cross- 
section views of an ultrasonic field is presented. The 
technique is applied to ultrasonic transducer radiation 
patterns with sample contour drawings of several cross- 
sections. A spatial frequency spectrxam formulation of 
diffraction theory is chosen for the computations which are 
performed with the aid of a fast Fourier transform routine. 

By experimentation an actual ultrasonic field was 
generated using a quartz transducer and the field was 
measured in a cross-sectional plane. A computer generated 
drawing of the measured data is compared to the drawing of 
a computer generated prediction of the field at the same 
plane . 

The computer model is applied to cin array of transducers 
to investigate the resulting ultrasonic field, as a pre- 
liminary use of the technique for array design. 
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INTRODUCTION 



I . 

The work in this study was in connection with ultrasonic 
imaging systems studies at the Naval Postgraduate School. 

The study undertaken here was generated by the need for a 
large-area ultrasonic field as delineated in Ref. 1. One 
possible way of producing a large planar wavefront is by 
synthesis of a large array from smaller transducers. Due to 
the large distribution of variables involved in designing a 
near-field array, a fast and practical means of producing 
near-field information was first required. A large body 
of literature is available for simple apertures in the near- 
field and C. M. McKinney ^t ^ recount the history of this 
research very concisely in Ref. 2. Exact expressions of the 
near-field are developed through essentially the same methods. 
They all involve finding a solution to the wave equation 
which satisfies certain assumed boundary conditions. The 
resulting integral equations are unwieldy even for the 
simplest aperture shapes because the field cannot be 
expressed in terms of elementary functions. The ability to 
project planes of data in the near-field to other planes 
within the near-field is therefore a most valuable facility 
in terms of making maximum use of both theoretical and 
empirical data available for design purposes. Developing a 
model to make these plane-to-plane projections was the major 
task undertaken in this study. A starting point for planning 
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the work was a study of the imaging technique reported by 
A. L. Boyer ^t _al [Ref. 3] . 

As a test of the model developed, the distribution of a 
square aperture's diffraction pattern was predicted at 
various planes throughout the field. The model was then 
used to predict the field of a quartz transducer. The 
actual field distribution of the quartz transducer modeled 
on the computer was measured after setting up the necessary 
experimental apparatus. A piezoelectric probe attached to 
a mechanical scan system was used to collect the data. A 
preliminary use of the data was the comparison of a predicted 
and a measured field distribution in a plane within the near- 
field of the source. The data taken in this experiment 
should be useful for future work in this and related areas 
of imaging system studies. 

Time allowed the use of the computer model to make pre- 
dictions of the near-field distribution of a transducer 
array. However, the use of the computer model for array 
design was done only in a preliminary manner by examining 
an assumed array which was relatively simple to model. This 
allowed the method of design to be evaluated although a 
design problem itself could not be accomplished within the 
time alloted for this work. 
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II. THEORY 



A. THE GENERAL SCHEME 

The purpose of this section is to present some of the 
fundamental concepts of an approach taken for computer 
modeling of the near-field of a transducer. It is not 
intended as a rigorous proof in the foixndations of spatial 
frequency spectrum theory [Ref. 4], but rather a justifica- 
tion of the computer scheme used in this study. The overall 
scheme is one which allows the prediction of a complex near- 
field distribution in a plane on the basis of a known complex 
field distribution in any other plane perpendicular to the 
direction of propagation. Based on the assumption that the 
planar distribution of a sound field at the transducer takes 
on the geometry and motion of the transducer itself, predic- 
tions of the sound field distribution at various planes 
throughout the near-field were made with the computer model. 
Figure 2-1 is a schematic which is helpful in organizing the 
operations involved in the following subsections. The sub- 
sections contain an explanation of the symbols which appear 
in Fig. 2-1. 

B. THE SPATIAL FREQUENCY SPECTRUM 

Let an isolated sound source be contained entirely within 
a source plane as shown below in Fig. 2-2. 
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Figure 2-1. The Propagation Model Sequence of Operations. 
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Figure 2-2. Two-dimensional Sound Source. 



Given that the source is a monochromatic pressure field in 
free space, and that x, y, and z are the cartesian coordi- 
nates of a typical field point, an expression of the pres- 
sure distribution in the source plane is 



P(x,y,z^t) = P(x,y,z^) g 

where t is time and w the angular frequency, 
are the following: 



( 2 . 1 ) 
Also defined 



P = the complex pressure phasor 
P = the peak pressure amplitude 

- the phase of the complex phasor 
z^ = the source plane position on the z axis 

For convenience an abbreviation in notion is used hereafter 
such that 



P^(x,y,t) = P(x,y,z,t) 



( 2 . 2 ) 
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where the subscript denotes the plane located at z. Assum- 
ing harmonic time dependence the expression for complex 
pressure is simplified to 

j4>„(x,y) 

® (2.3) 

In viewing P (x,y) as a two dimensional distribution for a 
z 

given z plane, the decomposition of the distribution into a 
spectrum of plane waves proves convenient. A decomposition 
of the distribution can be formed by use of a two-dimensional 
Fourier transform as given by Goodman [Ref. 5] . 

00 

P^(f^,fy) = J Jp^(x,y)exp[-j2Tr(f^x+f^y) ] dxdy (2.4) 

Goodman further shows that the original distribution P (x,y) 

z 

can now be written in terms of its spatial frequency spectrum 
by use of the inverse Fourier transform as shown by 

00 

^ 2 ^^'^"^ = fj ^z^^x'^y^ exp[ j2TT(f^x+f^y) ] df^df^ (2.5) 

The spatial frequencies, f^ and f^, can be related to a 
\mique orientation of the propagation vector, k, for a given 
plane wave in 3-dimensional space. Where k, the propagation 
constant, is given by k = 2it/X = 27rf/c, and f is the spatial 
frequency of a propagating plane wave in the direction of 
propagation; the vector, k, has magnitude, k, in the direc- 
tion of an outward normal from the planewave wavefront. Let 
the vector, ic, be separated into its components in the coor- 
dinate system by forming its dot product with a radius 
vector, r, such that 
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( 2 . 6 ) 



ic*r = k x+k y+k z = 2ii (f x+f y + f z) 

X y-^ z X y-' z 

Eq. 2.6 can be written in terms of the x-y plane spatial 
frequencies alone by making the following substitution: 




In the form given in Eq. 2.5, P (x,y) expresses the distri- 

z 

bution of sound pressure over the observed plane, z, as a 

superposition of plane waves and evanescent waves. Plane 

2 2 1 

waves are implied when f + f _< — « and evanescent waves are 

X y 

implied when f^ + f^ > [Ref. 6] . 

X y ^2 

C. PROPAGATION OF THE SPATIAL FREQUENCY SPECTRUM 

The propagating waveform must satisfy the homogeneous 

wave equation in 3-dimensional space. P (x,y) in the frequen- 

z 

cy spectrum formulation as given in Eq. 2.5 is required to 
satisfy the Helmholtz equation 

V^P + k^P = 0 (2.7) 

z z 



The solution of the Helmholtz equation leads to a second 
order differential equation in the variable z. Its solution 
is given in Ref. 5 as 



exp 



j2Tr(f x+f y) 

X y-' 



OO I • ’■■■■■ 

P„(x,y)= ffp'(f ,f )exp j2TT/-^ - f^-f^ 
z Id ^ ^ Y J ^ ^ 

( 2 . 8 ) 

Comparing Eq. 2.5 and Eq. 2.8 it is apparent that Eq. 2.8 
can be considered the two-dimensional Fourier transform of 
that portion of the expression under the integral ^ven by 



df df 
X y 
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) = exp j2TTZ 



y 




This expression suggests that propagation could be modeled 
in the Fourier transform domain using as a transfer function 



where z in Eq. 2.9 is taken as the distance between the 
reference plane and the observed plane. 

D. TWO-DIMENSIONAL LINEAR SYSTEM 

It remains to be shown that the propagation from one 
plane to another can be considered an overall linear system 
to justify the transfer function approach to modeling. In 
a linear system the total response is equal to the combina- 
tion of elementary responses, or mathematically stated the 
superposition property holds [Ref. 7]'. In order to demon- 
strate this property of frequercy spectrum propagation, it is 
convenient to express the waveform at an observed plane 
z = z - in terms of the waveform at a reference plane, z = z . 



To formulate a .wave distribution at the z = z plane as 

o 

shown in Fig. 2-3, a svibstitution of f into Eq. 2.8 and a 



regrouping of that expression leads to a planar distribution 
given by • • , 




(2.9) 



by merely taking the inverse transform of P^(f^,f^) • H(f^,f^) 



z 



<30 




( 2 . 10 ) 
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Figure 2-3. Boundaries and Distributions of a 
Propagated Wave Distribution at 
Two Planes. 



When the Fourier transform is denoted by F{ } and the 
inverse Fourier transform by F }, Eq. 2.10 becomes 



z o o 
o 



-1 



P ' (f f ) e 
z X y 
o 



j^^rf z 
z o 



( 2 . 11 ) 



Taking the Fourier transform of Eq. 2.11 and reordering the 
equation gives 



which can be equated to solve for P' (f ,f ) as given by 

^o ^ y 




At the plane z 



f ) = exp[-j27rf z ] ff P (x ,y ) 
y £' >• J z o J J z 0 -^0 

— oo O 

X exp [-j 2ir (f^x^+f^y^) ] dxdy (2.13) 

= z^ write the expression for P^ ^y 
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substituting Eq. 2.13 for P^(f^,f^) in Eq. 2.8. The result- 
ing expression for P terms of P (x,y) is: 

Z j Cl Cl z 

a o 

OO 00 

, <='aya> = / [ ®=Pt-3 2>tf2(2a-^o' // <Vo> 

Cl ..OO — OO O 



X exp[-j2TT (f^x^+f^y^) ] 



X exp[j2TT (f^x^+f^y^+f^z^) Idx^dy^df^df^ (2.14) 

By interchanging the order of integration of this equation 
and by making the simplification in notation that 



■^ao'^'a-^o'^a-yo* = // exp{j2z 



+ f (y^-y ) +f (z^-z ) ]} df df (2.15) 

y -^d -^o z d o J X y 



the final form of the expression of P at the z = z, plane 

Z ^ Cl 

d 

in terms of at the z = z^ plane can be expressed as 



^Z = //^Z 'Vo”'do'='d-=‘o'ya-i'o>'>Vya. 



Systems whose output can be expressed in this form as a 
function of the input are recognized to be linear systems 
[Ref. 8]. Thus the wave distribution at the z = z^ plane 
is given in terms of the distribution at the z = z^ plane by 
Eq. 2.16 when the system transfer function is given by 
Eq. 2.15. 

In the computer model developed for Eq. 2.9 evanescent 
waves were neglected and the transfer fvinction was therefore 
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bandlimited to those frequencies within a circle of radius 
1/X . This at first appears to limit the class of pressure 
distributions which can be considered to those which have 
bandlimited spatial frequency distributions. However, evan- 
escent waves decay exponentially with increased distance in 
the z-direction, and their effect on the field is insignifi- 
cant in .planes greater than a few wavelengths from the source. 
The region of the field where the evanescent waves have a 
significant effect on the pressure distribution is sometimes 
called the very-near-f ield. Acoustic imaging systems 
usually reconstruct the image from the pressure distribution. 
Reconstruction of an image requires that the evanescent 
waves be accounted for in the superposition integral of the 
type given in Eq. 2.16. Since the evanescent waves do not 
propagate, the near-field and far-field distributions are 
independent of their presence. Neglecting the evanescent 
waves, therefore, did not restrict the ability to make 
plane-to-plane predictions in the near-field as long as the 
reference distribution was measured at a distance greater 
than a few wavelengths from the transducer. Expressing 
Eq. 2.9 in a more appropriate form for the model which 
neglects the evanescent waves 



H(f ,f ) 

X y 



exp j27TZ -ir - 

Jx^ X y 




(2.17) 



0 



otherwise 
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E. THE DISCRETE FOURIER TRANSFORM 

Let X and Y be the sample intervals in the spatial x and 
y coordinate directions respectively. Let and be the 
corresponding sample intervals in the spatial frequency 
plane. The discrete Fourier transform (DFT) over a region - 
N and M samples long in the x and y directions respectively 
is defined as 



N-1 M-1 

F(kQ ) = E E f(nX,mY) 

^ ^ n=0 m=0 



where 

and 



X exp [- j (nXkfi^+mY£fi^) ] 



0 < k < N-1 and 0 < i < M-1 



X NX 



fj = 
y MY 



(2.18) 



These equations are in the form of those which were imple- 
mented with the fast Fourier transform (FFT) algorithm 
called "F0UR2" [Ref. 9]. The Fortran statements of these 
programs are included in Appendix A. The F0UR2 subroutine 
was used to perform the transformation to and from the 
spatial frequency spectrum domain. 



F. PROPAGATION BY TRANSFER FUNCTION 

When Eg. 2.17 is placed into a discrete format as 
follows : 

exp 






j2TTZ 



(kfi 

X v ^ z 



elsewhere 



(2.19) 



15 













4 







I 



I 



I 



I 



fi 











I 





or more simply written 



exp[j2irz ~ - if (k^+£-^) (Xfi ) ^ < 1 

■V / 1 ^ xy xy 



H 



k,i 



if (XQ ) ^ > 1 

xy — 



( 2 . 20 ) 



The propagation algorithm as given in Appendix A was used to 
perform the multiplication of the transformed input data by 
the transfer function (Eq. 2.20). 



G. ANALOG- TO-DIGITAL CONVERSION 

When predicted patterns of the waveform at various planes 
were calculated on the basis of their aperture function, the 
data was arranged in the shape of the aperture with uniform 
amplitude and phase samples everywhere within the aperture. 
However, the analog data collected in the experimental por- 
tion of this analysis had to be converted to sampled data 
before being utilized in the propagation program. The 
analog- to-digi tal , A/D, conversion process was done by the 
simultaneous utilization of an XDS 9300 digital computer and 
a COMCOR Ci 5000 analog computer. A lengthy description of 
the facilities and equipment used, the A/D conversion pro- 
cess, and the seven-to-nine track conversion was a necessary 
step in making the digital data of the XDS 9300 computer 
compatable with the IBM 360/67 computer requirements. 

H. AN INTERPOLATION OF DIGITAL DATA DISPLAY 

The subroutine CONTUR [Ref. 11] included as part of 
Appendix B, was used to produce separate analog displays of 
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the amplitude and phase at a desired plane. In this routine 
a pseudo three dimensional effect is presented by contour 
lines. Coarse analysis of an entire cross-section of the 
field can be performed from contour drawings. 

I. FLOWCHART OF STEPS IN COMPUTER OPERATIONS 

A flowchart of the computer operations in the order 
performed and a description of the equipment used in con- 
junction with each step is given in Fig. 2-4. 

J. INHERENT ADVANTAGES AND DISADVANTAGES OF THE DFT AND FFT 
Since Cooley and Tukey first reported the algorithm 

which they developed for the FFT [Ref. 12], the amount of 
literature which has been published on applications of this 
technique of computation and its advantages and disadvantages 
would fill volumes. Much insight into the FFT technique can 
be gained from Ref s» 13-15. Much less has been written 
about applying this technique for processing of data in the 
two dimensional spatial frequency domain. Reference 13 was 
the most germane discussion found for this study. It 
discusses the problems commonly called leakage and aliasing 
in articles which treat time- frequency analysis problems 
[Ref. 14] . 

The problem of leakage can sometimes work as an advan- 

,( 

tage in spatial frequency analysis. Leakage is an approxi- 
mation error which arises from a finite record length 
representation of an infinite length record. Its affect on 
the desired spectrum is equivalent to an undesired convolu- 
tion of a data window with the true frequency spectrum. The 
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Figure 2-4. Flowchart of Computer Operations for 
Propagation Modeling. 
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window occurs because of the geometry of the termination of 
the record. In this study a transducer array was in reality 
a truncated source and propagation of the data window 
resulted in a diffraction effect of the window geometry. 
Diffraction was modeled using this effect when the square 
shape of the sampled lattice was suitable. When the square 
window did not adequately model the geometry desired, the 
simplest method of modeling the edge diffraction was found 
to be one of leaving a border of zeroed data around the data 
field to be transformed. This method reduced the diffraction 
effects of the square sample lattice at the cost of a greatly 
increased memory requirement to model a given transducer. A 
discussion of the truncation approach to handling various 
geometries is given in Ref. 17. 

Aliasing is a problem associated with modeling a con- 
tinuous Fourier transform (CFT) by a discrete Fourier trans- 
form (DFT) . It is, therefore, an inherent disadvantage of 
the DFT. The term "aliasing" refers to high frequency 
components of the spatial function which are incorrectly 
represented due to widely spaced sampling. The effect is 
worse than the loss of the high frequency components since 
it results in erroneous coefficients in the lower fre- 
quencies of the spectrum. The effect can be avoided by 
sampling the waveform with the optimum spacing [Ref. 18] . 

The sample spacing when determined by the two-dimensional 
sample theorem frequently led to excessive numbers of 
samples necessary to represent a given transducer. Since 
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computations were computer memory limited, a corresponding 
limit was placed on the total number of samples available 
to represent a cross-section of the field. A compromise in 
choosing sample spacing and border width around tbe repre- 
sentative data was necessary. 

The disadvantages of the FFT associated with modeling 
the diffraction equations were insignificant when compared 
to those of the DFT. That is to say that the errors in 
modeling of the diffraction came from making the transition 
from a CFT to a DFT not in the FFT algorithm itself. The 
worth of the algorithm is measured in terms of speed, 
efficiency, memory required and round-off error. Of these, 
the requirement for a large memory is the only significant 
disadvantage. Although this study was seriously limited by 
the computer memory available, there are means to circumvent 
this problem by using more computation time. This is a 
matter of optimizing the algorithm to the particular region 
of the propagating waveform that is of interest. For this 
work, three versions of the algorithm listed in Appendix B 
were required. The subroutine "F0UR2" makes very effective 
use of time in performing the calculations of the FFT, how- 
ever, the total number of sample points (record length) must 
be a power of 2. Reference 19 gives an algorithm which 
can transform many different record lengths at a moderate 
cost in increased computation time. Implementation of this 
"mixed radix" form of algorithm is strongly recommended for 
future studies of transducer sound field analysis. 
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III. EXPERIMENTAL INVESTIGATION 



Experimental data was taken for two primary reasons. 
First was the desire for a means of evaluating the computer 
predictions. Second was the desire for actual data vice the 
assumed data at the face of the transducer as mentioned in 
Section II to evaluate the assumptions of the model itself. 

A mechanical scan probing system was constructed to measure 
the complex sound field of a transducer to be described. 

A. THE WATER TANK 

The water tank used in this experiment was constructed 
of plywood and lined with a one inch thickness of poly- 
styrene. The inner dimensions provided for a column of 
water 4 J 5 inches deep, 4 inches wide and 45 inches long. A 
cross section of this column was not large enough to con- 
sider the propagating wave entirely unbounded. However, it 
was found by experiment that the difference in nearby rigid 
and nearby release surface measurements were not significant 
until the hydrophone probe was within approximately ten 
wavelengths of the boundary. All data was taken with this 
consideration limiting the field of data. The length of the 
tank used was considered adequate to approximate an infinite 
medium at the frequency used. The tank was lined with poly- 
styrene which was a composition of pliable thermoplastic 
with a high absorption coefficient. Its impedance is very 
close to that of air for which p C = 415 rayls. 
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B. THE HYDROPHONE PROBE 



In design of the hydrophone (Fig. 3- 7a) , an attempt was 
made to minimize reflections affecting the free field 
measurement. The active face of the probe is 0.04 inches in 
diameter (0.83 wavelengths 0 1012 kHz). A thickness vibrator 
type receiver of PZT-5 ceramic composition is produced by 
Valpy-Fisher . .It is encapsuled with a bronze impedance 
matched backing and provided with a shielded electrical 
connector by the manufacturer. This type probe was mounted 
for use as a hydrophone in a plexiglass holder and coated 
with neopreme. 

C. THE TRANSDUCER 

The source transducer is shown in Fig. 3-lb. The crystal 
shown is a 2 inch by 2 inch square transducer of X-cut quartz. 
This transducer is designed as a thickness vibrator polished 
£>r a fundamental frequency of 1 MHz with good overtone opera- 
tion characteristics. The peak output signal from the probe 
was obtained when driving the transducer at 1012 kHz. Con- 
sequently that frequency was chosen for all data collected. 

The small offset from the transducer's fundamental resonance 
was necessary due to the mis~n'atch in resonant frequencies 
of transducer and probe. 

D. THE MECHANICAL SCAN EQUIPMENT 

The mechanical scan equipment moiinted above the measure- 
ment tank is shown in Fig. 3-2. The scanning system became 
cin important concern in fabricating the measurement apparatus 
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Figure 3-1. (a) The 2 Inch Square Quartz Crystal and Transducer Mounting, 

(b) Valpy-Fisher PZT-5 Probe and Hydrophone Mounting. 
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Figure 3-2. Measurement Tank and Mechanical Scan System. 





because three rigid system requirements depended upon the 
scan mechanism. The first was the requirement to insure 
that phase information be recorded accurately to within ±10 
degrees. This meant that the recording plane would have to 
be flat to within 0.04 millimeters at the frequency used of 
1012 kHz. A similar precaution was taken by Boyer et al as 
noted in Ref. 3. The second requirement was that the scan 
system be equipped to measure position accurately to within 
a small fraction of a wavelength. Still a third requirement 
was that it provide a means of taking data rapidly enough to 
collect a large quantity of data in a reasonable length of 
time. The heart of the scan system was a mechanical scan 
unit taken from a disassembled fire control computer. This 
unit, driven by a variable speed motor, positioned the probe 
horizontally across the width of the tank. A manually 
operated slide assembly positioned the probe in the vertical 
dimension to vary its depth. The coordinate position of the 
probe within the cross-sectional plane was calculated by the 
advance per turn of the threads on the positioning mechanisms. 
The number of turns was recorded by hand from a mechanical 
counter in the vertical, and by a potentiometer in the hori- 
zontal direction. Analog records were taken for each depth 
as the motor driven assembly scanned the width of the tank. 
Uniform depth increments were used between records so that 
ultimately a square matrix of equally spaced data points 
would result. The assembled system produced data which shows 
a phase variation across the face of the transducer of 20 
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degrees peak to peak. The phase information recorded in 
Fig. 3-3 was taken with the probe moving horizontally across 
the width of the transducer 9.1 wavelengths in front of it, 
at a point midway between the upper and lower edges of the 
transducer. There is an overall shift in the mean phase 
from the left edge to right edge of the transducer of 
approximately 30 ^degrees. This .indicates a misalignment of 
the transducer and recording planes by 0.12 millimeters 
maximum. Phase records were used in this manner to align 
the transducer to the recording plane in both the horizontal 
and vertical aspect. The maximum deviation in phase about 
the mean was ±10 degrees when taken directly in front of the 
transducer. On the basis of this and similar records it was 
estimated that phase measurements were accurate to within 
±10 degrees at the frequency used of 1012 kHz. Turns of the 
scan positioning mechanisms were recorded accurately to 
within 1/36 of a turn with a resultant position accuracy of 
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Figure 3-3. Phase Variation across the Face of a Quartz 
Transducer Radiating at 1012 kHz. 
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approximately 1/50 of a wavelength. The sample spacing was 
9/10 of a wavelength at uniformly spaced intervals. 

E. ELECTRONIC MEASURING EQUIPMENT 

The general layout of electronic measurement equipment 
is shown in Fig. 3-4 and a photograph of the physical layout 
is given in Fig. 3-5. The electronic equipment layout was 
an assembly of standard pieces of equipment. A 
general description and specif ication summary is given 
in Appendix B for the major components. 

F. EXPERIMENTAL PROCEDURE 

It was pertinent to the data reduction process that the 
diameter of the hydrophone face was less than the sample 
spacing, thus avoiding overlapping of the samples. Had 
imaging been the desired use of the data collected, it would 
have been ideal to sample the field at half-wavelength 
intervals. The scanning system described in this experiment 
could have been used for imaging if a more complex data 
reduction process were used [Ref. 3] . The distance between 
the transducer and the recording plane was adjusted by 
moving the transducer which was mounted on a traveling 
holder. Separation of transducer and probe was measured by 
a scale graduated in 1/64 inch. If time had allowed repeated 
experiments of this type, a mechanically driven threaded 
shaft and traveler system would have been used for convenience 
in observing ultrasonic wave behavior in the direction of 
propagation. 
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Figure 3-4. Electronic Equipment Layout. 
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Figure 3-5. Electronic Equipirent Setup. 



G. TYPICAL EXPERIMENTAL RESULTS 



Figure 3-6 shows six typical records as recorded at 
successive depths by the analog type scanning equipment. 
Amplitude, phase, and potentiometer outputs are recorded 
simultaneously on separate channels. The experiment was 
designed to generate data for a 64 by 64 sample lattice, so 
the potentiometer output has a spike corresponding to each 
of 64 equally spaced positions of the hydrophone along the 
width of the record. This channel was used in the A/D 
conversion to trigger the sample circuit. Figure 3-7 
shows the analog playback of the sampled data as recorded 
on magnetic tape after the A/D conversion. Comparison of 
the measured record and the converted record shows good 
agreement. 
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Figure 3-6. Six Typical Records of Amplitude, Phase, and Potentiometer Output 
(reading top to bottom) as Recorded at Successive Depths by the 
Scanning Equipment. 
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Figure 3-7. Six Typical Records of the Sampled Data Records Showing Amplitude 
and Phase (reading top to bottom) as Played Back from a Magnetic 
Tape in Analog Presentation. 



IV. EXPERIMENTAL vs THEORETICAL RESULTS 



A test of the computer program was the generation of a 
far-field diffraction pattern. This was done by using a 
small aperture whose diffraction pattern could be considered 
to be in the far-field at distances well within the near- 
field of the sample lattice used. The effect of the sample 
lattice was introduced in the preceding section in the 
discussion of the DFT where the scheme of making the ^erture 
to be modeled much smaller than the "data window" was pro- 
posed. This scheme allows increased latitude in the range 
of z for which the diffraction model can be considered valid. 
If a sufficient number of discrete samples are taken to 
accurately represent a square aperture, and a large border 
of zeroed data is left around the aperture, then the window 
does not disturb the desired diffraction pattern. It was 
convenient to consider this effect in the following light: 

The very- near- fie Id of an aperture is essentially a 
parallel projection of the aperture. The diffraction of a 
small aperture diverges very rapidly and the far-field is 
reached rather quickly in moving in the z-direction. The 
diffraction pattern of a large aperture diverges more slowly 
and the far-field is reached at much greater distances of 
propagation. (The same radiation frequency for both aper- 
tures is assumed.) The terminology "near-field of the data 
window" will be used to refer to the region in which the 
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modeled aperture is valid. (Equation 2.9 which is the exact 
and continuous form of the propagation model has no such 
limitation on the region of the z axis in which it is valid. 
The limitation is a result of a finite sample lattice and 
with a large enough sample space the diffraction of any 
aperture can be projected into the far-field with this 
model. ) 



W. K. Pratt and H. C. Andrews [Ref. 20] have generated 
computer far-field diffraction patterns which were used for 
comparison to test the computer program for this work. 

Figure 4-1 shows a typical diffraction pattern of a single 
square aperture. The amplitude variations of the field 
(ranging from 0.0 to 3.1 units) are present in five contour 
levels, for a plane 96 wavelengths (96X) from the aperture. 
The divisions between levels are 0.7 units apart. The 
width of the aperture, a, was 6X and the width of the data 
field was 38X across. The pattern shown v;as as expected 
from the analytic solution for a square aperture [Ref. 5]. 

A variety of criteria are found in the literature for deter- 
mining the transition from the near-field to the far-field 
region. According to Goodman [Ref. 5], this transition can 
be given as 



z » 



2 2 
k (x“^+y^) 



max 



(4.1) 



where z and k are as previously defined and x and y are 
measured from the center of the aperture as shown in figs. 
2.2 and 2.3. Writing Eq. 4.1 in terms of sampled data for 
a square lattice gives the criterion 
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Figure 4-1. 



Far-field Diffraction Pattern (amplitude) of a 
Square Aperture. a/A = 6 , z = 96X. 
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(4.2) 



z » 




X 



From this equation the maximum z for which the window does 
not appreciably interfere with the observed diffraction 
pattern is taken as 



5tt .^2. 
z = N X 
max 12 



(4.3) 



For example, the pattern shown in Fig. 4.1 is a square 
lattice of 32 samples giving a far-field for the data window 
beginning at z 1340X. Applying the same equation (Eq. 4.1) 
to the maximum x and y, where x = y = N/16, the far- 
field of the aperture is where 

2tt ,N„\2 

z » (ipT- X) or z >> 17. 4X 

A 16 

It was expected, therefore, that Fig. 4-1 would be a reason- 
ably accurate presentation of the far-field pattern of a 
square aperture out to z = 1340X. 

In the preceding discussion the word aperture was used 
instead of the word transducer because of the type of distri- 
bution assumed at the transducer plane. When the following 
assumptions are made the transducer diffraction pattern is 
based entirely on its aperture function: 

1. The motion of the transducer has uniform amplitude 
over its entire face. 

2 . The motion of the transducer is in phase over the 
entire face. 

t 

3. The motion of the medium immediately adjacent to the 
face of the transducer has the same motion as the transducer 



itself. 
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4. The medium outside of the edges of the transducer 
does not move due to the motion of the transducer. 

The program was tested by its generation of aperture diffrac- 
tion patterns. The ability to predict transducer diffraction 
patterns on the basis of an aperture diffraction pattern will 
be discussed in light of the experimental data. 

The transducer shown in Fig. 3-lb -was modeled on the 
IBM 360 computer and the field at z = 73X was predicted 
based upon its square aperture function. Figures 4-2 and 
4-3 are the predictions of amplitude and phase respectively. 
Note that the transducer was a 2 inch square sheet of quartz, 
while the active face assumed in the model was a 1 3/4 inches 
square aperture. This was a correction for that portion of 
the transducer obscured by the mount. Additionally, the 
restriction in motion offered by the mount at the edges 
violates the assumption of uniform amplitude. This effect 
can be seen from the amplitude record made at 9. IX in front 
of the transducer. Figure 4-4 shows a typical amplitude 
measurement at z = 73 X. The assumed shape at the trans- 

ducer is shown by the heavy dotted line superimposed over 
the data record. Figures 4-5 and 4-6 are contour plots of 
the sampled data at z = 9. IX as recorded after the A/D con- 
version. The disparities in the predicted and measured 
contour patterns is attributable to the following inaccuracies: 

1. The wavefront and measurement planes were slightly 
misaligned as seen from the phase record of Fig. 3-3. In 
order for the fourfold symmetry of Figs. 4-2 and 4-3 to be 
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Figure 4-2. Predicted Near-field Diffraction Pattern 
(amplitude) of a Square Transducer. 
a/X = 30, z = 73X. 
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Figure 4-3. Predicted Near-field Diffraction Pattern 
(phase) of a Square Transducer. 
a/A = 30, z = 73A. 
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Figure 4-4. Amplitude Variation Across the Face of a Quartz 
Transducer Radiating at 1012 kHz for a/X = 30, 




Figure 4-5. Measured Near-field Diffraction Pattern 

(amplitude) of a Square Transducer a/X = 30, 
z = 73X. 
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Figure 4-6. Measured Near-field Diffraction Pattern (phase) 
of a Square Transducer a/X = 30, z = 73X. 



present in the measured pattern the source must be aligned 
symmetrically v/ith respect to the measurement plane. The 
misalignment also caused the propagating waveform to shift 
its center with respect to the center of the sample lattice. 

2. Data were not taken of the entire plane. To con- 
serve time in the recording and processing of data, only half 
of the plane was sampled; symmetry was assumed and the half- 
plane data was folded about a horizontal through the center. 
This gave a false symmetry to the data in view of the align- 
ment difficulty. 

3. Recording errors in the measurements taken were on 
the order of ±10 degrees for phase and ±4% for amplitude, 
due to inaccuracies of the recording equipment. The A/D 
conversion process introduced additional error, but no 
information was available on the accuracy of the program 
used other than the visual comparison from displays such 
as the one shown in Fig. 3-7. 

4. Inaccuracies in the free- field assumption caused 
by the probe hydrophone, and by the boundaries in the tank 
were neglected. 

5. Inaccuracies were inherent in evaluating a trans- 
ducer's diffraction pattern based on its aperture function 
alone. 

It is impossible to evaluate the validity of modeling 
the transducer by its aperture function on the basis of the 
comparisons presented here. If the data recorded at 9. IX 
were propagated to the plane at 73X and compared to the 
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measurements made there, the composite effect of the inaccu- 
racies listed in (1) to (4) above could be determined. It 
was concluded for the purpose of this experiment that the 
aperture function model was sufficient for modeling small 
transducers, and was found to be an adequate model when the 
ratio of a/A was less than 10. 

Many acoustic imaging systems of interest require planar 
ultrasonic fields at frequencies of 1 MHz and higher. An 
array design approach to supplying this need was the motiva- 
tion for developing the computer prediction capability. As 
a preliminary study to the design problem a 3 x 3 array of 
square transducers was chosen. After some experimenting with 
the program it was found that at least 25 samples per trans- 
ducer was required to produce an acceptable far-field 
diffraction pattern for a square aperture. Each square 
transducer was represented by an aperture 6 wavelengths on 
a side. The spacing between these transducers had to be 
wide enough so that interference between adjacent transducers 
did not occur until a distance where the far-field pattern 
of the individual transducers had already developed. This 
effect is evident in the pattern shown in Fig. 4-7, where 
the sound field can be described as a sum of nine spatially 
separated far-field patterns. When the array was propagated 
to greater and greater distances, the far-field patterns 
expanded and merged together in a way which made the indi- 
vidual far-field pattern distinguishable in the overall 
field pattern. When the array was propagated to that distance 
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Figure 4-7. Diffraction Pattern for a transducer Array of 
9 Uniformly Spaced Elements. The elements are 
squares, 6 wavelengths on a side. The pattern 
is a prediction for a plane 60 wavelengths 
from the source plane. The contour levels are 
2.8 units apart. 
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Figure 4-8. A Predicted Diffraction Pattern for a Transducer 
Array of 9 Uniformly Spaced Elements. The 
figure is one quadrant of a pattern which has 
four-fold symmetry. The elements are squares 
6 wavelengths on a side. The pattern resents a 
plane 600 wavelengths from the source plane. 

The contour levels are 0.28 units apart. 
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where adjacent elements began to interfere, the field 
diffraction pattern became very complex and all similarity 
to the far-field pattern of a square was lost. Instead 
the pattern looked somewhat like the near-field diffraction 
pattern of a finely sampled square aperture. Figure 4-8 
is an example of this effect which the author has called 
"the’^near-f ield diffraction of -far-field sources". This 
nomenclature alludes to the fact that the individual trans- 
ducers were in the far-field according to the criteria of 
Eq. 4.1 while the outer dimensions of the composite array 
placed the source as a whole in the near-field. This figure 
is the result of contouring only one quadrant and using a 
128 X 128 sample lattice. (Note the symmetry of Fig. 4-8, 
and the improved resolution gained by the larger number of 
samples used.) 

The use of the computer model appears to have much 
promise for designing an array to produce a desired waveform. 
There are a large number of variables which determine the 
form of the array diffraction. It was found, however, that 
the number of trials which must be run can be significantly 
narrowed by defining dimensionless qiantities which model 
many different combinations of element size, element spacing 
and element shape. Variables which were not looked into in 
this preliminary study were the number of elements in the 
array, and the possible variations in driving the individual 
elements in various phase and amplitude combinations. If 
all feasible combinations of the variables mentioned were 
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to be examined for a design, it would be exhaustively time 
consuming and expensive. A major advantage of this approach 
is the intuitive feeling the user gains by running a few 
simple cases using well known aperture shapes. Familiarity 
with the geometry of some two-dimensional Fourier trans- 
forms [Ref. 21] coupled with a feeling for the general trends 
of propagation in the spatial frequency domain, would allow 
the designer to eliminate a large class of possible variables 
and to concentrate on those which have some intuitive promise. 
Facilities are available at the Naval Postgraduate School 
for immediate graphical response to a change in variables 
using the graphics display terminal which presents data on 
a rectangular cathode-ray tube. This type of data presenta- 
tion used with a routine similar to the CONTUR subroutine 
would make the large number of variables involved manageable. 
The IBM/360 version of the propagation program as presented 
in Appendix B requires 250K core storage with the CONTUR 
subroutine included. This version requires approximately 
56 seconds of model 67 cpu time to compile and calculate 
the data for the first projection of a square matrix of 
64 X 64 samples. Each additional propagation of the data 
.to another plane takes approximately 20 seconds. The time- 
memory requirement of the program appears to be adaptable to 
use with the graphic display units. 
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V. CONCLUSIONS 



A recently developed technique for imaging by reproduc- 
tion from object diffraction patterns has been applied to 
analyzing the sound source itself. The technique is not 
limited to ultrasound, but high frequency imaging applica- 
tions were of particular interest for this study. Conse- 
quently, the experimental measurements were conducted at 
1 MHz. Although most of the concepts and arguments can be 
applied to frequencies well below IMHz, such applications 
were not considered herein. 

A technique of computer modeling of the diffraction 
equations has been emphasized. The capability of the 
computer facilities at the Naval Postgraduate School to 
perform such diffraction calculations in a fast concise 
format has been accented. 

The use of this technique in designing an array was a 
driving incentive in the work reported, and from the results 
of the preliminary array model examined it appears to be a 
worthwhile approach to production of a planar field of 
ultrasound. 

The experimental portion of this work produced data for 
am elementary comparison with a computer prediction. The 
data has many applications for future work to give infor- 
mation about the transducer used, about the capability of 
the computer program, and about imaging system capabilities. 
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APPENDIX A 



PROPAGATION PROGRAM - IBM 360/67 COMPUTER 

THIS PROGRAM MODELS PROPAGATION OF AN ULTRASONIC SOUND 
SOURCE AND PERFORMS PL AN E-TO-PLAN E PROJECTIONS OF A TWO- 
DIMENSIONAL COMPLEX PRESSURE DISTRIBUTION. THE REFERENCE 
PLANE DATA IS ENTERED INTO TWO REAL ARRAYS OF 64 X 64 
SAMPLES GIVING A TOTAL RECORD LENGTH OF 4096 POINTS. THE 
TWO REAL ARRAYS ARE USED TO FORM ONE COMPLEX ARRAY. THE 
COMPLEX ARRAY IS TRANSFORMED TO THE SPATIAL FREQUENCY 
PLANE USING THE SUBROUTINE F0UR2. PROPAGATION TO ANOTHER 
PLANE IS PERFORMED BY A TRANSFER FUNCTION ALGORITHM WHICH 
IS A MULTIPLIER IN THE TRANSFORMATION DOMAIN. THE INVERSE 
TRANSFORM IS TAKEN USING THE SUBROUTINE F0UR2 AND THE DATA 
AT THE NEW PLANE IS DISPLAYED ON THE PAGE PRINTER AND ON 
THE CALCOMP PLOTTER BY USE OF THE SUBROUTINE CONTUR. 

THIS IS ONE OF THREE BASIC PROGRAMS WRITTEN TO OPTIMIZE 
THE TIME-MEMORY TRADE-OFF FOR A PARTICULAR TRANSDUCER TYPE. 

SEE *A COMPUTER TECHNIQUE FOR NEAR-FIELD ANALYSIS OF AN 

ULTRASONIC TR AN SDUCE R* ,GRI GG S , C. A. , M. S . THESIS, U.S. NAVAL 
POSTGRADUATE SCHOOL, CAL I FORNIA, MONTEREY, 1973 



IMPLICIT C0MPLEX*8(D) 

REAL*8 TITLE(12) 

REAL*8 LABEL(12) 

DISTANCES TO BE LABELED ON CONTOUR PLOTS ARE ENTERED HERE 

REAL+8 ZTITLE(l)/' OX •/ 

LOGICAL=f:l LTG(3) 

DIMENSION CLM (4) , CLP (4) 

DIMENSION DATA(64,64J , NN(2 ) ,DATAVC(4096 ) , AMAGVC(4096) 
1PHAVC( 4096), XVEC( 64) , YVEC{ 64) , YVECP(64) 
C0MM0N/AREA1/AMAG(64,64) ,PHA(64, 64) 

DIMENSION DC0PY(64,64) ,DC0PVC(4096 ) 

EQUIVALENCE (DATA( 1, 1) , DATA VC ( 1 ) ) , ( AMAG( 1 , 1 ) , AMAGVC(l) 
EQUIVALENCE ( PHA ( 1 , 1 ) , PHAV C ( 1 J ) 

C 

C INITIAL VALUES OF N AND M 

DATA NN/64,64/ 

LABELS WHICH APPEAR ON THE CONTOUR DRAWINGS ARE ENTERED 
HERE WITH DISTANCES DEPENDENT UPON DATA ABOVE AND NOT ON 
THE ACTUAL VALUE OF 'Z* IN THE PROGRAM 

DATA TITLE/' • , 

1 • A* , 

2 'MPLITUDE', 

3 • CONTOUR', 

4 

5 ' ', 

5 'ONE SQ T', 

6 'RANSDUCE', 

7 'R Z=', 

9 * 

C ' C , 

C 'A GRIGGS'/ 

DATA LABEL/' ', 

1 • ' , 

2 ' PHASE', 

3 * CONTOUR', 

4 

5 ' ' , 

5 'ONE SQ T', 

6 'RANSDUCE', 

7 »R Z=', 

C ' ', 

C • C , 

C 'A GRIGGS'/ 
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LTG(1)=.TRUE. 

LTG(2 }=. FALSE. 

LTG( 3) =.FALSE. 

NINT=NN(1) 

MINT=NN(2 ) 

NREL=(NINT+2) /2 
MREL=(MINT+2)/2 
REALN=FLOAT(NI NT) 

REALM=FLOAT(MINT) 

NDIM=NN(1) 

MDIM=NN(2) 

NVEC = NDIM=^'MDIM 
B=0.0 

PI=3. 14159265 
DJ=( 0.0, 1.0) 

DO 999 K=1,NVEC 
DATA VC (K)=(0. 0,0.0) 

AMAGVC(K)=(0.0) 

PHAVC(K)=(0.0) 

DCOPVC(K) =(0.0) 

999 CONTINUE 
C 

C READING MEASURED DATA 
C 

AS AN OPTIONAL METHOD OF ENTERING DATA SIMULATED DATA CAN 
BE GENERATED AND ENTERED DIRECTLY INTO THE COMPLEX ARRAY. 

READ (5,100) (AMAGVC(K),K=1,2048) 

READ(5,100) (PHAVC(K) ,K=1,2048) 

100 FORMAT (5(F9. 4, 5X) ) 

C 

C READING DATA INTO DATA(I,J) 

C 

DO 1001 J=l,29 
DO 1001 1=1,64 
PEAK=(AMAG(I,J)+8.0)/9.6 
PHI = PHA( I , J)='‘( PI/99 .0 ) 

CLIP=J+3 

DATA( I, CL IP)=CMPLX(PEAK, PHI ) 

1001 DCOPY( I,CLIP)=OATA(I,CLIP) 

C 

C CNE RECORD TAKEN OUT OF ORDER 
C 

THIS PARTICULAR VERSION PERFORMED SOME CORRECTIONS TO THE 
DATA BECAUSE OF THE DATA RECORDING PROCESS. TWO ROWS IN 
THE MATRIX WERE INTERCHANGED AND THE DATA WAS SHIFTED TO 
CORRECT FOR THE POINT OF SYMMETRY. IT ALSO FOLDS A HALF 
PLANE OF DATA ASSUMING TWO-FOLD SYMMETRY. 

J = 28 

DO 1002 1=1, 64 
0RDER=J+1 

DATA (I , J) =DCOPY( I , ORDER) 

1002 CATA( I ,ORDER)=DCOP Y( I , J) 

DO 1003 J=33,64 

DO 1003 1=1,64 
F0LD=65-J 

1003 CATA( I , J)=DATA( I, FOLD) ‘ 

CALL CONVRT(DATAVC,AMAGVC, PHAVC) 

c : 

WRITE(6,3130) 

3130 FORMATCl* ,T57,' INPUT ARRAY*/) 

CALL OUTPUT 
C 

SINCE THE DISCRETE TRANSFORM IS INDEXED FROM ZERO TO N MINUS 
ONE THE TRANSFORM IS CENTERED ABOUT A POINT SHIFTED BY 
N DIVIDED BY TWO SAMPLES TO THE RIGHT AND UP IN THE FIRST 
QUADRANT THE INPUT DATA MUST BE SHIFTED ALSO TO MATCH 
THIS AMOUNTS TO A PHASE SHIFT IN THE SPATIAL FREQUENCY PLANE 
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DO 2001 I=1,NDIM 
DO 2001 J=1,MDIM 
DEX=(0J=<'PI=^( I + J) J 
DSHIF=CEXP(DEX) 

DATA( I , J) = DATA(I,J J’S'DSHIF 

2001 CONTINUE 
C 

CALL F0UR2IDATA.NN, 2,-1) 

DATA! I, J) ARE NOW THE FOURIER COEFFICIENTS (SHIFTED) 

DO 2002 K0PY=1,NVEC 

2002 DCOP VC { KOPY) =DATAVC( KOPY) 

Z=0.0 

L00P=1 

ALAM=1. 08853774 
DO 200 3 KAT=1,NVEC 

2003 DATAVC(KAT)=DCOPVC{KAT) 



C PROPAGATION BY SPATIAL FREQUENCY TRANSFER FUNCTION. 



WAVE = 1.0/( ALAM4ALAM) 
FREQX=1.0/REALN 
FREQY=1.0/REALM 
DO 3003 K=1,NDIM 
.KINT=K-(NINT+2) /2 
PRIMEK=FLOAT(KINT) 
FXK=FREGX* PRIM5K 
FXKSQ=FXK*FXK 
DO 3003 L=1,MDIM 
LINT=L-(MINT+2) /2 
PRIMEL=FLOAT(L INT) 

FYL = FREQY=^^PRIMEL 

FYLSQ = FYL=f-FYL 

EVANS* (WAVE-FXKSQ-FYLSQ) 

3000 RADI AN=2 .0=<=PI 
ROOT=SQRT( EVANS) 
DPR0P=DJ:5:RADIAN*Z*R00T 
DTRANS=CEXP(DPROP) 

DATA(K,L )=DATA( K,L)=5'DTRANS 
GO TO 3002 

3001 DATA(K,L)=(0. 0,0.0) 

3002 CONTINUE 

3003 CONTINUE 



CALL F0UR2 ( DATA,NN, 2, 1) 

DO 6004 I=1,NDIM 
DO 6004 J=1,MDIM 
PDIM = 1.0/{ REALN=«'REALM) 
DIMPLX=CMPLX(POIM,B) 
DEX=(-CJ*PI*( I+J) ) 
DSHIF=CEXP(DEX> 

DATA{ I, J )=DATA( I , J )*DIMPLX*DSHIF 
6004 CONTINUE 

CALL CONVRT(DATAVC ,AMAGVC, PHAVC) 
WRITE(6,6005) Z,REF 
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6005 FORMAT ( • 1 • , T5 7 , • INPUT ARRAY PROPAGATED TO Z=*,F5.1/ 
2*0* ,T57, ’PHASE REFERENCE =’,F5.1,' DEGREES’/) 

CALL OUTPUT 

C REFERENCE PHASE TO UNDIFFRACTED RAY IN DEGREES 

DO 6009 K=l, NVEC 

PHAVC(K) = ( PH AVC (K)* (180, 0/PI )) +180.0 
REF=PHA(NREL,MREL) 

PHAVC (K)=ABS( PHAVC(K )-REF) 

IF(PHAVC(K).LE. 180.0) GO TO 6009 
PHAVC (K)=360.0-P HA VC(K) 

6009 CONTINUE 

TITLE( 10)=ZTITLE(L00P) 

LABEL( 10)=ZTITLE(L00P) 

CALL CCNTUR ( AMAG , 64 ,64,64 , CLM, -4t T ITL E, 9, 9 , LTG) 

CALL CONTUR ( PHA, 64, 64, 64, CLP , -4 ,LABEL ,9 ,9 , LTG) 

Z=Z+20 .0 

STOP 

END 



SUBROUT INE CONVRT ( DAT AVC,AMAGVC, PHAVC) 

IMPLICIT C0MPLEX=*'8 (D) 

DIMENSION DATAVC(4096) ,AMAGVC( 4096) ,PHAVC(4096) 
DO 1080 K=l,4096 
REALP=REAL (DATAVC( K) ) 

AIMAGP=AIMAG(DATAVC( K) ) 

AMAGVC (K)=SQRT ( REALP-^2+AIMAGP**2) 

TRY=AMAGVC (K) 

IF(TRY.EQ.O.O) GO TO 1070 
PHAVC ( K) =ATAN2 (AIMAGP, REALP) 

GO TO 1080 
1070 PHAVC(K)=0.0 
1080 CONTINUE 
RETURN 
END 



C SUBROUTINE OUTPUT FOR A 64 X 64 



SUBROUTINE OUTPUT 

C0MM0N/AREA1/AMAG(64,64) ,PHA(64, 64) 
WRITE( 6,3031) (J,J=33,64) 

3031 FORMAT( ’0’,4X,32I4) 

DO 3040 1=33,64 

WRITE! 6,3032) ( I ,(AMAG{ I , J) , J=33 ,64) ) 

3032 FORMAT ( ’O’ , 13, 1X,32F4. 1) 

WRITE (6, 3033) { PHA( I,J) ,J=33,64) 

3033 FORMAT! • O’ ,4X,32F4.1/) 

3040 CONTINUE 

RETURN 

END 
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SUBROUT INE FOUR2( DATA , NN , ND I M, I SIGNJ 
DIMENSION DATA(1),NN(1) 

IF(NDIM-l) 700, 1, 1 

1 NT0T=2 

DO 2 IDIM=1,NDIM 
IF(NN( IDIM))700,700,2 

2 NTOT=NTOT*NN( lOIM) 

TW0PI=6. 2831 853070 
RTHLF=. 7071067812 
NP1=2 

DO 600 IDIM=1,NDIM 
N=NN{ I DIM) 

NP2=NP1=^N 
IF(N-1 )700,600,100 
100 NP2HF=NP2/2 
J = 1 

DO 160 I2=1,NP2,NP1 
IF( J-I2) 110, 130, 130 
110 IlMAX=I2+NPl-2 

DO 120 11 = 12, UMAX, 2 
DO 120 13=11, NT0T,NP2 
J3=J+I3-I2 
TEMPR=DATA(I3) 

TEMPI=DATA( 13+1) 

DATA! I3)=DATA( J3 ) 

DATA! 13+1) =0ATA( J3+1) 

CATA( J3)=TEMPR 
120 CATA(J3+1)=TEMPI 
130 M=NP2HF 

140 IF(J-M ) 160, 160, 150 
150 J=J-M 
M=M/2 

IF(M-NPl) 160, 140, 140 
160 J=J+M 

NP1TW=NP1+NP1 

IPAR=N 

310 IF( I PAR-2) 350,330,320 
320 IPAR=IPAR/4 
GO TO 310 

330 DO 340 11=1, NPl, 2 

DO 340 K1=I1,NT0T,NP1TW 

K2=K1+NP1 

TEMPR=DATA(K2) 

,TEMPI=DATA(K2+1) 
DATA(K2)=DATA(K1)-TEMPR 
DATA(K2+1)=DATA( Kl + D-TEMPI 
DATA(K1)=DATA(K1)+TEMPR 
340 DATA(K1+1 )=DATA(K1+1)+TEMPI 
350 MMAX=NP1 

360 IF(MMAX-NP2HF) 370, 600, 600 
370 LMAX=MAX0(NPlTW,MMAX/2) 

DO 570 L=NP1, LMAX, NPITW 
M=L 

IF(MMAX-NP1)420,420,380 
380 THETA=-TWOPI=?'FLOAT(M) /FLOAT! 4*MMAX) 
IF ( IS IGN)400, 390, 390 
39C THETA=-THETA 
400 WR=COS(THETA) 

WI=SIN(TH£TA) 

410 W2R=WR*WR-WI*WI 
W2I=2.0*WR^WI 
W3R = W2R^'WR-W2I*WI 
W3I = W2R*WI + vg2I*WR 
420 DO 530 11=1, NPl, 2 
KMIN=I PAR*M+I1 
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IF (MMAX-NP 1)430,430,440 
430 KMIN=I1 
440 KDIF=I PAR*MMAX 
450 KSTEP=4*KDIF 

DO 520 K1=KMIN,NT0T,KSTEP 

K2=K1+KDIF 

K3=K2+KDIF 

K4=K3+KDIF 

IF(MMAX-NP1)460,460,480 
460 U1R=DATA(K1)+DATA( K2) 

U1 I=DATA(K1+1 )+DAT A(K2+li 
U2R=DATA{K31+DATA( K4 ) 
U2I=DATA(K3+1)+DATA(K4+1) 
U3R=CATA(K1)-DATA(K2) 
U3I=DATA(K1+1)-DATA( K2+li 
IF(ISIGN)470,475,475 
470 U4R=0ATA(K3+1)-DAT A(K4 + 1) 
U4I=0ATA( K4)-DATA( K3) 

GO TO 510 

475 U4R=DATA{ K4+1)-DATA(K3+1) 
U4I=DATA(K3)-0ATA(K4) 

GO TO 510 

480 T2R=W2R=i'DATA(K2)-W2I’«'DATA( K2+1 ) 
T2I=W2R=<=DATA(K2+1)+W2I*DATA(K2) 
T3R=WR*DATA(K3 )-W I*DATA(K3 + 1) 
T3I = WRVDATA( K3 + 1 ) + WI^i'DATAI K3) 
T4R=W3R*DATA(K4)-W3I=^DATA( K4+1) 
T4 1 = W3 R^DATA ( K4+1 ) +W3 1 =^‘DAT A ( K4 ) 
U1R=DATA(K1)+T2R 
U1I=DATA(K1+11+T2I 
U2R=T3R+T4R 
U2I=T3I+T4I 
U3R=DATA(K1)-T2R 
U3I=DATA(K1+1)-T2I 
IF( ISIGN)490,500,500 
490 U4R=T3I-T4I 
U4I=T4R-T3R 
GO TO 510 
500 U4R=T4I-T3I 
U4I=T3R-T4R 
510 DATA(K1)=U1R+U2R 
DATA(K1+1)=U1I+U2I 
DA TA(K2)=U3R+U4R 
CATA(K2+1)=U3I+U4I 
DATA (K3 1=U1R-U2R 
DATA(K3+1)=U1I-U2I 
DATA(K4)=U3R-U4R 
520 DATA(K4+1) =U3I-U4I 
KMIN=4=^(KMIN-I 1)+I 1 
KDIF=KSTEP 

IF(KDIF-NP2HF)450,450,530 
530 CONTINUE 
M=M+LMAX 

IF(M-MMAX)540,540,570 
540 IF( ISIGN) 550,560,560 
550 TE.MPR=WR 

WR=( WR+WI)*RTHLF 
WI=(WI-TEMPR)>^RTHLF 
GO TO 410 
560 TEMPR=WR 

WR= (WR-Wn*RTHLF 
WI=(TEMPR+WI 1*RTHLF 
GO TO 410 
570 CONTINUE 

IPAR=3-IPAR 
MMAX=MMAX+MMAX 
GO TO 360 
600 NP1=NP2 
700 RETURN 
END 
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APPENDIX B 



GENERAL DESCRIPTION AND SPECIFICATION SUIIMARY 
OF MAJOR ELECTRONIC COMPONENTS 



EQUIPMENT 

NAME 


TAPE 

RECORDER 


POWER 

AMPLIFIER 


VECTOR 

VOLTMETER 


MANUFACTURER 
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Instruments 

for 
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Hewlett 

Packard 
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8405A 


GENERAL 
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D.C. recorder 
output has 
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VoFo/^e occiput 
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ancC pko-se- 


INPUT / OUTPUT 
CHARACTER- 
ISTICS 


In^t level : o.s 4o 
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VTlodxcLcLt-( 
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+ lO dbrvi stayidl- 
(XViL . 

I n p ot 1 m pedai^ee 
S~0 obrns 
Oa+pot trhpcJAMcc- 
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"Boj^pass tipple 
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Yollojoe ; 

Ck. A.- I.s mV +o 
1 V 

CkBt lOO^V +t> 
1 V rrv\s 
Ou4 Oo*t 

Ar%!ipli4aJL^: O 'Lo 
ivd_.a. i 4% 

: O i o 

± o.s VdL.c. t 6% 



56 



LIST OF REFERENCES 



1. Powers, J. P., and Griggs, C., "Large-Area Ultrasonic 

Fields for Acoustic Imaging," Ultrasonics Symposium 
Proceedings , p. 104-108, October 1972. 

2. McKinney, C. M. , Harvel, K. W. , and Ellis, G. E., 

"Characteristics of Line and Disk Underwater Sound 
Transducers in the Near-to Farfield Transition 
Region," The Journal of the Acoustical Society of 
America , v. 51, p. 1076-1081, 17 August 1971. 

3. Boyer, A. L. , and others, "Teconstruction of Ultrasonic 

Images by Backward Propagation," Acoustical Holo- 
graphy , V. 3, p. 333-348, 31 July 1970. 

4. Ratcliffe, J. A., "Some Aspects of Diffraction Theory 

and Their Application to the Ionosphere," Reports 
on Progress in Physics , v. 19, The Physical Society, 
London, 1956. 

5. Goodman, J. W. , Introduction to Fourier Optics , p. 48-50, 

62-65, McGraw-Hill, 1968. 

6. Shewell, J. R. , and Wolf, E., "Inverse Diffraction and a 

New Reciprocity Theorem," Journal of the Optical 
Society of America, v. 58, p. 1596-1603, December 
1968. 

7. Cooper, G. R. , and McGilem, C. D., Methods of Signal 

and System Analysis , p. 11-12, Holt, Rinehart, and 
Winston, Inc., 1967. 

8 . Thomas , J . B . , An Introduction to Statistical Communica - 

tion Theory , p. 140-144, Wiley, 1969. 

9. Brenner, N. M. , "Three Fortran Programs that Perform the 

Cooley-Tukey Fourier Transform," IEEE Transactions 
on Audio and Electroacoustics, v. AU-17, No. 2, 
p. 93-103, June 1969. 

10. Jones, R. D. , Time Series Analysis of Analog Data by 

Analog-to-Digital and Digital Processing Methods at 
the Naval Postgraduate School , M.S. Thesis, Naval 
Postgraduate School, California, Monterey, 1971. 

11. Naval Postgraduate School Technical Note No. 0211-03, 

Plotting Package for NPGS IBM 360/67 , by P. C. 
Johnson, February 1969. 



57 



12. Cooley, J. W. , and Tukey, J. VI., "An Algorithm for the 

Machine Calculation cf Complex Fourier Series," 
Mathematics of Computers, v. 19, p. 297-301, April 
1965. 

13. Sondhi, M. M. , "Reconstruction of Objects from Their 

Sound-Diffraction Patterns," The Journal of the 
Acoustical Society of America, v. 46, no. 5, 
p. 1158-1164, 30 April 1969. 

14. Bergland, G. D. , "A Guided Tour of the Fast Fourier 

Transform," IEEE Spectrum , p. 41-52, July 1969. 

15. Cochran, W. T. and others,’ "What is the Fast Fourier 

Transform?", p. 45-55, IEEE Transactions on Audio 
and Electroacoustics , v. AU-15, no. 2, June 1969. 

16. Bertram, S., "Frequency Analysis Using Discrete Fourier 

Transform, " IEEE Transactoins on Audio and Electro - 
acoustics, V. Au-18, no. 4, p. 495-500, December 
1970. 

17. Peterson, D. P. and Middleton, D. , Information Control, 

V. 5, p. 279-323, 1962. 

18. Papoulis, A., Systems and Transforms with- Applications 

in Optics , p. 119-132, McGraw-Hill, 1968. 

19. Singleton, R. C. , "An Algorithm for Computing the Mixed 

Radix Fast Fourier Transform," IEEE Transactions 
on Audio and Electroacoustics , v. AU-17, no. 2, 
p. 93-103, June 1969. 

20. Pratt, W. K. , and Andrews, H. C. , "Computer Calculated 

Diffraction Patterns," Applied Optics , v. 7, no. 2, 
p. 378-379, February 1968. 

21. Bracewell, R. , The Fourier Transform and Its Applications , 

p. 246-247, McGraw-Hill, 1965. 



58 



INITIAL DISTRIBUTION LIST 



No. 



1. Defense Documentation Center 
Cameron Station 
Alexandria, Virginia 22314 

2. Library, Code 0212 
Naval Postgraduate School 
Monterey, California 93940 

3. Asst. Professor J. P. Powers, Code 52 Po 
Department of Electrical Engineering 
Naval Postgraduate School 

Monterey, California 93940 

4. Asst. Professor A, I. Eller, Code 61Er 
Department of Physics 

Naval Postgraduate School 
Monterey, California 93940 

5. LCDR Carlton A. Griggs, USN 
USS WAHOO (SS 565) 

FPO San Francisco, California 96601 

6. Asst. Professor V. M. Powers, Code 52Pw 
Department of Electrical Engineering 
Naval Postgraduate School 

Monterey, California 93940 



Copies 

2 

2 

1 

1 

1 

1 



59 





I 



UNCLASSIFIED 

Secuntv Classification 



DOCUMENT CONTROL DATA • R & D 

fSfcurity classification of titio, body of abstract and indexing annotation must be entered tvhen f/te overall report Is cfasslfled) 



^CinaTinC activity (Corporate author) 

Naval Postgraduate School 
Monterey, California 93940 



2«. REPORT SECURITY CLASSIFICATION 

Unclassified 



2b. CROUP 



■Report title 

A COMPUTER TECHNIQUE FOR NEAR-FIELD ANALYSIS 
OF AN ULTRASONIC TRANSDUCER 



descriptive notes (Type of report and^incluaive dates) 

Master's Thesis; March, 1973 



rAUTMORiS) (First namSt middle initial, last name) 



Carlton A. Griggs; Lieutenant Commander, United States Navy 



HEPORT DATE 

March 1973 



ll. CON TRAC T OR GRANT NO. 



PROJECT NO. 



7a. total no. of paces 



61 



7b. NO. OF REFS 



21 



9a. originator*! report NUMBER(S) 



9b. OTHER REPORT NO(S) (Any other numbera that may be aeelgc^ed 
this report) 



10. OISTRI BU TION STATEMENT 



Approved for public release; distribution unlimited. 



M. 5UPPLEMEN T A RV NOTES 



12. SPONSORING MILITARY ACTIVITY 



Naval Postgraduate School 
Monterey, California 93940 



13. ASSTRAC T 



A computer technique for obtaining near-field cross- 
section views of an ultrasonic field is presented. The 
technique is applied to ultrasonic transducer radiation 
patterns with sample contour drawings of several cross- 
sections. A spatial frequency spectrum formulation of 
diffraction theory is chosen for the computations which 
are performed with the aid of a fast Fourier transform 
routine. 

By experimentation an actual ultrasonic field was 
generated using a quartz transducer and the field was 
measured in a cross-sectional plane. A computer generated 
drawing of the measured data is compared to the drawing 
of a computer generated prediction of the field at the 
same plane. 

The computer model is applied to an array of trans- 
ducers to investigate the resulting ultrasonic field, 
as a preliminary use of the technique for array design. 



DD ; 



'.r..i473 

S/N 0l0l*807>6ai I 



(PAGE I) 



60 



UNCLASSIFIED 



Security CUttification 



A- 3140 S 



UNCLASSIFIED 

Security Clatsification 



KEY WORDS 


LINK A 


LINKS 


LINK C 


ROLE 


W T 


ROLE 


W T 


ROLE 


W T 


ULTRASONICS 
ACOUSTIC IMAGING 
TRANSDUCER DESIGN 

COMPUTER DIFFRACTION CALCULATION 





























DD .'r‘..1473 

$/H 0I0I-«07'«82I 



UNCLASSIFIED 



61 



Security Clataification 



A> 3 t 409 



14353S 



Thes i s 

G829 

c.l 



Thes i s 

G 829 Griggs 
c.l A computer technique 

for near-field analysis 
of an ultrasonic trans- 
ducer. 



Griggs 



143536 



A computer technique 
for near-field analysis 
of an ultrasonic trans- 
ducer. 



thesG829 

A computer technique for near-field anal 




3 2768 001 03801 1 
DUDLEY KNOX LIBRARY 




